{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018076d6-8f21-4949-82cc-52ea73f366d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from linearmodels.panel import PanelOLS\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import LabelEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2c54d73-7e8f-4c96-810c-d914e82b3520",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"\\df_1.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06406552-7e5c-4707-ac84-5ee4cc5dd2eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a02ebe92-283a-4983-a1e9-e7a2ec67f03f",
   "metadata": {},
   "source": [
    "# 1. Analysis 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d457bb0-f9f1-47fe-a496-7a020b2d856e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming df is your DataFrame and it's already loaded\n",
    "df = df.set_index(['si_gun_gu', 'year'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af76d6cd-b03b-4dd9-a35a-f3b3dbe22057",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb37d830-049c-49df-ab3e-4c6c803ea83c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    X = sm.add_constant(X)\n",
    "    model = PanelOLS(y, X, entity_effects=True)\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    return results\n",
    "\n",
    "# Define independent variables\n",
    "independent_vars = ['Urbanization', 'Income', 'Homeownership', 'Age', 'Sex']\n",
    "\n",
    "# Models for different dependent variables\n",
    "results_happy_evi = run_fe_regression(df, 'Happy', ['EVI'] + independent_vars)\n",
    "results_happy_ndvi = run_fe_regression(df, 'Happy', ['NDVI'] + independent_vars)\n",
    "results_relaxed_evi = run_fe_regression(df, 'Relaxed', ['EVI'] + independent_vars)\n",
    "results_relaxed_ndvi = run_fe_regression(df, 'Relaxed', ['NDVI'] + independent_vars)\n",
    "results_energetic_evi = run_fe_regression(df, 'Energetic', ['EVI'] + independent_vars)\n",
    "results_energetic_ndvi = run_fe_regression(df, 'Energetic', ['NDVI'] + independent_vars)\n",
    "results_anxiety_evi = run_fe_regression(df, 'Anxious', ['EVI'] + independent_vars)\n",
    "results_anxiety_ndvi = run_fe_regression(df, 'Anxious', ['NDVI'] + independent_vars)\n",
    "results_sad_evi = run_fe_regression(df, 'Sad', ['EVI'] + independent_vars)\n",
    "results_sad_ndvi = run_fe_regression(df, 'Sad', ['NDVI'] + independent_vars)\n",
    "results_depressed_evi = run_fe_regression(df, 'Depressed', ['EVI'] + independent_vars)\n",
    "results_depressed_ndvi = run_fe_regression(df, 'Depressed', ['NDVI'] + independent_vars)\n",
    "results_angry_evi = run_fe_regression(df, 'Angry', ['EVI'] + independent_vars)\n",
    "results_angry_ndvi = run_fe_regression(df, 'Angry', ['NDVI'] + independent_vars)\n",
    "results_stressed_evi = run_fe_regression(df, 'Stressed', ['EVI'] + independent_vars)\n",
    "results_stressed_ndvi = run_fe_regression(df, 'Stressed', ['NDVI'] + independent_vars)\n",
    "results_tired_evi = run_fe_regression(df, 'Tired', ['EVI'] + independent_vars)\n",
    "results_tired_ndvi = run_fe_regression(df, 'Tired', ['NDVI'] + independent_vars)\n",
    "results_lonely_evi = run_fe_regression(df, 'Lonely', ['EVI'] + independent_vars)\n",
    "results_lonely_ndvi = run_fe_regression(df, 'Lonely', ['NDVI'] + independent_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c6c6bb4-6223-43cc-a33c-bf43d7f73a86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Extract for each model\n",
    "variables = ['EVI', 'NDVI']\n",
    "results_list = [results_happy_evi, results_happy_ndvi, results_relaxed_evi, results_relaxed_ndvi,\n",
    "                results_energetic_evi, results_energetic_ndvi, results_anxiety_evi, results_anxiety_ndvi,\n",
    "                results_sad_evi, results_sad_ndvi, results_depressed_evi, results_depressed_ndvi, results_angry_evi,\n",
    "                results_angry_ndvi, results_stressed_evi, results_stressed_ndvi, results_tired_evi, \n",
    "                results_tired_ndvi, results_lonely_evi, results_lonely_ndvi\n",
    "]\n",
    "\n",
    "data_for_plot = []\n",
    "for var in variables:\n",
    "    for i, result in enumerate(results_list):\n",
    "        if i % 2 == (0 if var == 'EVI' else 1):  # To ensure we match models with corresponding variables\n",
    "            coefs, ci_lower, ci_upper = extract_coefs(result, var)\n",
    "            data_for_plot.append([var, f'Model {(i//2)+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3620240-b2d7-4807-ab0d-64e57f73c02d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(15, 8))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.2\n",
    "positions = {'EVI': -offset, 'NDVI': 0}\n",
    "colors = {'EVI': 'green', 'NDVI': 'blue'}\n",
    "\n",
    "# Plotting each variable\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    plt.errorbar(adjusted_positions, subset['Coef'], yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                 fmt='o', label=var, capsize=5, markersize=13, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        plt.text(adjusted_positions[j], subset['Coef'].iloc[j] + 1.5, f'{subset[\"Coef\"].iloc[j]:.3f}', \n",
    "                 ha='center', va='bottom', fontsize=16, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "plt.axhline(0, color='red', linewidth=3, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Model 1\\n(DV: Happy)', 'Model 2\\n(DV: Relaxed)', 'Model 3\\n(DV: Energetic)', \n",
    "                 'Model 4\\n(DV: Anxious)', 'Model 5\\n(DV: Sad)', 'Model 6\\n(DV: Depressed)',\n",
    "                 'Model 7\\n(DV: Angry)', 'Model 8\\n(DV: Stressed)', 'Model 9\\n(DV: Tired)', \n",
    "                 'Model 10\\n(DV: Lonely)']\n",
    "plt.xticks(range(len(custom_labels)), custom_labels, size=15)\n",
    "plt.yticks(size=16)\n",
    "\n",
    "# Increasing legend text size\n",
    "plt.legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large')\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab58e731-8c39-4963-b91e-b4a990c00be2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77aaded0-baa1-4ecb-a344-f6e5d76d7202",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    X = sm.add_constant(X)\n",
    "    model = PanelOLS(y, X, entity_effects=True)\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    return results\n",
    "\n",
    "# Define independent variables\n",
    "independent_vars = ['Urbanization', 'Income', 'Homeownership', 'Age', 'Sex']\n",
    "\n",
    "# Models for different dependent variables\n",
    "results_positive_evi = run_fe_regression(df, 'Positive_feelings', ['EVI'] + independent_vars)\n",
    "results_positive_ndvi = run_fe_regression(df, 'Positive_feelings', ['NDVI'] + independent_vars)\n",
    "results_negative_evi = run_fe_regression(df, 'Negative_feelings', ['EVI'] + independent_vars)\n",
    "results_negative_ndvi = run_fe_regression(df, 'Negative_feelings', ['NDVI'] + independent_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9575db04-0749-4ab9-a357-badfb959c3db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Extract for each model\n",
    "variables = ['EVI', 'NDVI']\n",
    "results_list = [\n",
    "    results_positive_evi, results_positive_ndvi, results_negative_evi,\n",
    "    results_negative_ndvi,\n",
    "]\n",
    "\n",
    "data_for_plot = []\n",
    "for var in variables:\n",
    "    for i, result in enumerate(results_list):\n",
    "        if i % 2 == (0 if var == 'EVI' else 1):  # To ensure we match models with corresponding variables\n",
    "            coefs, ci_lower, ci_upper = extract_coefs(result, var)\n",
    "            data_for_plot.append([var, f'Model {(i//2)+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91869dde-42a6-4e7d-a798-407556300b0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(9, 5))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.2\n",
    "positions = {'EVI': -offset, 'NDVI': 0}\n",
    "colors = {'EVI': 'green', 'NDVI': 'blue'}\n",
    "\n",
    "# Plotting each variable\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    plt.errorbar(adjusted_positions, subset['Coef'], yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                 fmt='o', label=var, capsize=5, markersize=17, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        plt.text(adjusted_positions[j], subset['Coef'].iloc[j] + 1.4, f'{subset[\"Coef\"].iloc[j]:.3f}', \n",
    "                 ha='center', va='bottom', fontsize=18, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "plt.axhline(0, color='red', linewidth=4, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Model 1\\n(DV: Positive Affectivity)', 'Model 2\\n(DV: Negative Affectivity)']\n",
    "plt.xticks(range(len(custom_labels)), custom_labels, size=16)\n",
    "plt.yticks(size=16)\n",
    "\n",
    "# Increasing legend text size\n",
    "plt.legend(title='Variable', fontsize='x-large', title_fontsize='x-large')\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d72d1df7-bca3-4815-899c-612928a842d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    X = sm.add_constant(X)\n",
    "    model = PanelOLS(y, X, entity_effects=True)\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    return results\n",
    "\n",
    "# Define independent variables\n",
    "independent_vars = ['Urbanization', 'Income', 'Homeownership', 'Age', 'Sex']\n",
    "\n",
    "# Define dependent variables\n",
    "dependent_vars = ['Positive_feelings', 'Negative_feelings']\n",
    "\n",
    "# Generate results for all combinations of dependent variables and green indices (EVI, NDVI)\n",
    "results = {}\n",
    "\n",
    "for dep_var in dependent_vars:\n",
    "    for green_var in ['EVI', 'NDVI']:\n",
    "        results[f'{dep_var.lower()}_{green_var.lower()}'] = run_fe_regression(df, dep_var, [green_var] + independent_vars)\n",
    "\n",
    "# Function to generate predictions\n",
    "def predict_effects(results, var_name, var_range, independent_vars):\n",
    "    predictions = []\n",
    "    for value in var_range:\n",
    "        data = {var_name: value, 'const': 1}\n",
    "        for var in independent_vars:\n",
    "            data[var] = df[var].mean()  # Using mean of other independent vars\n",
    "        X_new = pd.DataFrame(data, index=[0] * len(df.index))\n",
    "        X_new.index = df.index  # Set the same index as the original data\n",
    "        X_new = X_new[['const'] + [var_name] + independent_vars]\n",
    "        pred = results.predict(exog=X_new)\n",
    "        predictions.append(pred.mean().item())  # Use mean to aggregate the predictions\n",
    "    return predictions\n",
    "\n",
    "# Ranges for predictions\n",
    "evi_range = np.linspace(0.1, 0.4, 100)\n",
    "ndvi_range = np.linspace(0.2, 0.7, 100)\n",
    "\n",
    "# Plotting the results\n",
    "fig, axes = plt.subplots(1, 2, figsize=(18, 8))\n",
    "\n",
    "titles = [\"(a) Positive Affectivity\", \"(b) Negative Affectivity\"]\n",
    "\n",
    "for i, dep_var in enumerate(dependent_vars):\n",
    "    ax = axes[i]\n",
    "    \n",
    "    predictions_evi = predict_effects(results[f'{dep_var.lower()}_evi'], 'EVI', evi_range, independent_vars)\n",
    "    predictions_ndvi = predict_effects(results[f'{dep_var.lower()}_ndvi'], 'NDVI', ndvi_range, independent_vars)\n",
    "\n",
    "    # Plot for EVI\n",
    "    ax.plot(evi_range, predictions_evi, label='EVI', color='green', linestyle='-', linewidth=4)\n",
    "    ax.text(evi_range[0], predictions_evi[0], f'{predictions_evi[0]:.2f}', fontsize=26, color='green', ha='right', va='bottom')\n",
    "    ax.text(evi_range[-1], predictions_evi[-1], f'{predictions_evi[-1]:.2f}', fontsize=26, color='green', ha='left', va='top')\n",
    "\n",
    "    # Plot for NDVI\n",
    "    ax.plot(ndvi_range, predictions_ndvi, label='NDVI', color='blue', linestyle='--', linewidth=4)\n",
    "    ax.text(ndvi_range[0], predictions_ndvi[0], f'{predictions_ndvi[0]:.2f}', fontsize=26, color='blue', ha='right', va='bottom')\n",
    "    ax.text(ndvi_range[-1], predictions_ndvi[-1], f'{predictions_ndvi[-1]:.2f}', fontsize=26, color='blue', ha='left', va='top')\n",
    "\n",
    "    ax.set_xlabel('Vegetation Index', fontsize=22)\n",
    "    ax.set_ylabel('Predicted Value', fontsize=22)\n",
    "    ax.legend(fontsize=25)\n",
    "    ax.set_title(titles[i], fontsize=30, pad=17)\n",
    "    ax.tick_params(axis='x', labelsize=18)\n",
    "    ax.tick_params(axis='y', labelsize=18)\n",
    "\n",
    "# Adjust layout for better spacing\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83e4a3f7-f215-4d38-bd61-a13aa20559dd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "226e4809-3dff-4dc3-a461-f8c333b74efa",
   "metadata": {},
   "source": [
    "## 2. Analysis 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23c6d734-36d0-4c58-af9a-0bd8e0bef9f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"\\df_2.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa88ee64-4e4c-46b2-9777-e30a4318c5e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55a4a072-38b5-485f-8a45-0d45733bbc3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "\n",
    "# Drop rows with missing values for simplicity\n",
    "df = df.dropna(subset=['Happy', 'Relaxed', 'Anxious', 'Sad', 'Depressed', 'Angry', 'Stressed', 'Tired',\n",
    "                       'Energetic', 'Lonely', 'EVI', 'NDVI', 'sex', 'marital', 'age', 'education', 'income_2', \n",
    "                       'homeownership', 'religion', 'year'])\n",
    "\n",
    "# Ensure proper encoding of categorical variables\n",
    "df['sido_1'] = LabelEncoder().fit_transform(df['sido'])\n",
    "df['si_gun_gu_1'] = LabelEncoder().fit_transform(df['si_gun_gu'])\n",
    "\n",
    "# Verify the grouping variable\n",
    "assert df['si_gun_gu_1'].isnull().sum() == 0, \"Grouping variable has missing values.\"\n",
    "assert len(df['si_gun_gu_1'].unique()) > 1, \"Grouping variable should have more than one unique value.\"\n",
    "\n",
    "# Function to run mixed effects model\n",
    "def run_mixed_model(data, formula, re_formula):\n",
    "    model = smf.mixedlm(formula, data, groups=data[\"si_gun_gu_1\"], re_formula=re_formula)\n",
    "    result = model.fit()\n",
    "    return result\n",
    "\n",
    "# Define dependent variables\n",
    "dependent_vars = ['Happy', 'Relaxed', 'Energetic', 'Anxious', 'Sad', 'Depressed', 'Angry', 'Stressed', 'Tired', 'Lonely']\n",
    "\n",
    "# Generate formulas\n",
    "formulas_evi = [f'{dv} ~ EVI + sex + marital + age + education + income_2 + homeownership + religion + C(year)' for dv in dependent_vars]\n",
    "formulas_ndvi = [f'{dv} ~ NDVI + sex + marital + age + education + income_2 + homeownership + religion + C(year)' for dv in dependent_vars]\n",
    "\n",
    "# Run mixed models\n",
    "results_evi = [run_mixed_model(df, formula, '~EVI') for formula in formulas_evi]\n",
    "results_ndvi = [run_mixed_model(df, formula, '~NDVI') for formula in formulas_ndvi]\n",
    "\n",
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Extract for each model\n",
    "variables = ['EVI', 'NDVI']\n",
    "data_for_plot = []\n",
    "for i, dep_var in enumerate(dependent_vars):\n",
    "    coefs_evi, ci_lower_evi, ci_upper_evi = extract_coefs(results_evi[i], 'EVI')\n",
    "    data_for_plot.append(['EVI', dep_var, coefs_evi, ci_lower_evi, ci_upper_evi])\n",
    "    \n",
    "    coefs_ndvi, ci_lower_ndvi, ci_upper_ndvi = extract_coefs(results_ndvi[i], 'NDVI')\n",
    "    data_for_plot.append(['NDVI', dep_var, coefs_ndvi, ci_lower_ndvi, ci_upper_ndvi])\n",
    "\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Dependent Variable', 'Coef', 'CI Lower', 'CI Upper'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c03ad9e0-03f9-4803-97d0-df78f79daae0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the results\n",
    "fig, axes = plt.subplots(2, 5, figsize=(12, 4))\n",
    "\n",
    "for i, dep_var in enumerate(dependent_vars):\n",
    "    ax = axes[i // 5, i % 5]\n",
    "    \n",
    "    subset_evi = df_plot[(df_plot['Variable'] == 'EVI') & (df_plot['Dependent Variable'] == dep_var)]\n",
    "    subset_ndvi = df_plot[(df_plot['Variable'] == 'NDVI') & (df_plot['Dependent Variable'] == dep_var)]\n",
    "    \n",
    "    # Plot for EVI\n",
    "    ax.errorbar(1, subset_evi['Coef'], yerr=[subset_evi['Coef'] - subset_evi['CI Lower'], subset_evi['CI Upper'] - subset_evi['Coef']],\n",
    "                fmt='o', label='EVI', capsize=5, markersize=10, color='blue')\n",
    "    ax.text(1, subset_evi['Coef'] + 0.05, f'{subset_evi[\"Coef\"].values[0]:.3f}', ha='center', va='bottom', fontsize=10, color='blue')\n",
    "\n",
    "    # Plot for NDVI\n",
    "    ax.errorbar(2, subset_ndvi['Coef'], yerr=[subset_ndvi['Coef'] - subset_ndvi['CI Lower'], subset_ndvi['CI Upper'] - subset_ndvi['Coef']],\n",
    "                fmt='o', label='NDVI', capsize=5, markersize=10, color='green')\n",
    "    ax.text(2, subset_ndvi['Coef'] + 0.05, f'{subset_ndvi[\"Coef\"].values[0]:.3f}', ha='center', va='bottom', fontsize=10, color='green')\n",
    "\n",
    "    ax.set_xticks([1, 2])\n",
    "    ax.set_xticklabels(['EVI', 'NDVI'], fontsize=12)\n",
    "    ax.set_ylabel(dep_var, fontsize=14)\n",
    "    ax.axhline(0, color='red', linewidth=1, linestyle='--')\n",
    "    ax.legend(fontsize=12)\n",
    "    ax.tick_params(axis='y', labelsize=12)\n",
    "\n",
    "# Adjust layout for better spacing\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "388dc61e-9790-4a08-a34d-2c6a3c2233f2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39ab04af-0399-49fc-8e72-79f99347af2b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
